
*-----------------------------------------------------------------------------
// Figure 4: The effects of road toll on car ownership, before and after the BEV expansion
*-----------------------------------------------------------------------------

*------------------------------
// Step 1: Run regressions and save estimates
*------------------------------

** Select variables to loop over: 
local outcomelist BEV_fam_yes car_fam_count car_fam_yes ICE_fam_count   
local treatvar1list toll_fam_mean       
local treatvar2list ptl_fam_km_mean    

*-------------------
/// Pre (2008-2010) 
*------------------

// loop over couples and singles
forvalues i=0(1)1 {

use    "${dataout}MainDataset" , clear
keep if year>=2008 & year<=2010
keep if couple==`i'

	if `i'==0 {
		//* Single-adult: Generate one value per row values 
		egen employed = rowmax(employed1 employed2) 
		egen retired = rowmax(retired1 retired2) 
		egen age = rowmax(age1 age2) 
		egen wies_decile = rowmax(wies1_decile wies2_decile) 
		egen wealth_decile = rowmax(wealth1_decile wealth2_decile) 
		egen grputd = rowmax(grputd1 grputd2) 
		egen time = rowmax(time1 time2) 
		egen dist = rowmax(dist1 dist2) 
		egen grk_bed_num  = rowmax(grk_bed1_num  grk_bed2_num) 
		egen PublicTransitTime  = rowmax(PublicTransitTime_1  PublicTransitTime_2) 
		egen VENTETID  = rowmax(VENTETID_1  VENTETID_2) 
		egen BOARDINGS  = rowmax(BOARDINGS_1  BOARDINGS_2) 
		egen TAKST  = rowmax(TAKST_1  TAKST_2 )
		
		local household     $householdS
		local employment    $employmentS 
		local age 		    $ageS
		local income 	    $incomeS 
		local education	    $educationS
		local distance 	    $distanceS   
		local time 		    $timeS  
		local FE		    $FES 
		local publictime	$publictimeS   
		local publicquality $publicqualityS
		local clustervar	$clustervarS
	}
	
	if `i'==1 {
		
		local household     $household
		local employment    $employment
		local age 		    $age
		local income 	    $income 
		local education	    $education
		local distance 	    $distance   
		local time 		    $time 
		local FE		    $FE
		local publictime	$publictime   
		local publicquality $publicquality
		local clustervar	$clustervar
	}


foreach outcome in `outcomelist'{
foreach treatvar1 in `treatvar1list'{
foreach treatvar2 in `treatvar2list'{

reghdfe `outcome' `treatvar1'  `treatvar2'  `age'    `distance'  `time'  `publictime'  `publicquality'   , ///
absorb(`FE'   `household'  `income'  `employment'  `education'  i.year   ) vce(cluster `clustervar') 

summarize `outcome' if e(sample)==1
estadd scalar MeanDep = r(mean), replace
summarize `treatvar1' if e(sample)==1
estadd scalar Treat1 = r(mean), replace
summarize `treatvar2' if e(sample)==1
estadd scalar Treat2 = r(mean), replace

estadd local year="2008-10" , replace
estadd local grkFE="\checkmark", replace
estadd local grkbFE="\checkmark", replace
estadd local HHcon="\checkmark", replace
estadd local HHWork="\checkmark", replace
estadd local PublicTime="\checkmark", replace
estadd local grkFE_W="Yes", replace
estadd local grkbFE_W="Yes", replace
estadd local HHcon_W="Yes", replace
estadd local HHWork_W="Yes", replace
estadd local PublicTime_W="Yes", replace

estimates save "${ster}reg_g4_2008_2010_O_`outcome'_T1_`treatvar1'_T2_`treatvar2'_couple`i'" , replace
}
}
}
}


*-------------------
/// Post (2015-2017) 
*------------------

// Re-use estimates created for Table 1 and Table 2

// loop over couples and singles
forvalues i=0(1)1 {

use    "${dataout}MainDataset" , clear
keep if year>=2015 & year<=2017
keep if couple==`i'

	if `i'==0 {
		//* Single-adult: Generate one value per row values 
		egen employed = rowmax(employed1 employed2) 
		egen retired = rowmax(retired1 retired2) 
		egen age = rowmax(age1 age2) 
		egen wies_decile = rowmax(wies1_decile wies2_decile) 
		egen wealth_decile = rowmax(wealth1_decile wealth2_decile) 
		egen grputd = rowmax(grputd1 grputd2) 
		egen time = rowmax(time1 time2) 
		egen dist = rowmax(dist1 dist2) 
		egen grk_bed_num  = rowmax(grk_bed1_num  grk_bed2_num) 
		egen PublicTransitTime  = rowmax(PublicTransitTime_1  PublicTransitTime_2) 
		egen VENTETID  = rowmax(VENTETID_1  VENTETID_2) 
		egen BOARDINGS  = rowmax(BOARDINGS_1  BOARDINGS_2) 
		egen TAKST  = rowmax(TAKST_1  TAKST_2 )
		
		local household     $householdS
		local employment    $employmentS 
		local age 		    $ageS
		local income 	    $incomeS 
		local education	    $educationS
		local distance 	    $distanceS   
		local time 		    $timeS  
		local FE		    $FES 
		local publictime	$publictimeS   
		local publicquality $publicqualityS
		local clustervar	$clustervarS
	}
	
	if `i'==1 {
		
		local household     $household
		local employment    $employment
		local age 		    $age
		local income 	    $income 
		local education	    $education
		local distance 	    $distance   
		local time 		    $time 
		local FE		    $FE
		local publictime	$publictime   
		local publicquality $publicquality
		local clustervar	$clustervar
	}



foreach outcome in `outcomelist'{
foreach treatvar1 in `treatvar1list'{
foreach treatvar2 in `treatvar2list'{

reghdfe `outcome' `treatvar1'  `treatvar2'  `age'    `distance'  `time'  `publictime'  `publicquality'   , ///
absorb(`FE'   `household'  `income'  `employment'  `education'  i.year   ) vce(cluster `clustervar') 

summarize `outcome' if e(sample)==1
estadd scalar MeanDep = r(mean), replace
summarize `treatvar1' if e(sample)==1
estadd scalar Treat1 = r(mean), replace
summarize `treatvar2' if e(sample)==1
estadd scalar Treat2 = r(mean), replace

estadd local year="2008-10" , replace
estadd local grkFE="\checkmark", replace
estadd local grkbFE="\checkmark", replace
estadd local HHcon="\checkmark", replace
estadd local HHWork="\checkmark", replace
estadd local PublicTime="\checkmark", replace
estadd local grkFE_W="Yes", replace
estadd local grkbFE_W="Yes", replace
estadd local HHcon_W="Yes", replace
estadd local HHWork_W="Yes", replace
estadd local PublicTime_W="Yes", replace

estimates save "${ster}reg_g4_2015_2017_O_`outcome'_T1_`treatvar1'_T2_`treatvar2'_couple`i'" , replace
}
}
}
}


*------------------------------
// Step 2: Load all estimates
*------------------------------

capt eststo drop *
capt estimates drop *

** Select variables to loop over: 
local yearlist  	  2015_2017  2008_2010  
local outcomelist     BEV_fam_yes  ICE_fam_count  car_fam_count car_fam_yes
local treatvar1list   toll_fam_mean     
local treatvar2list   ptl_fam_km_mean    

foreach outcome in `outcomelist'{
foreach treatvar1 in `treatvar1list'{
foreach treatvar2 in `treatvar2list'{
foreach y in `yearlist' {

estimates use "${ster}reg_g4_`y'_O_`outcome'_T1_`treatvar1'_T2_`treatvar2'_couple1" 
eststo c_`y'_`outcome'

estimates use "${ster}reg_g4_`y'_O_`outcome'_T1_`treatvar1'_T2_`treatvar2'_couple0" 
eststo s_`y'_`outcome'

}
}
}
}

estimates dir


*------------------------------
// Step 3: Plot coefficients
*------------------------------

local treatvarlist   toll_fam_mean   
local hhtypelist c s

foreach hhtype in `hhtypelist'{
foreach treatvar in `treatvarlist'{
	
	
if "`treatvar'"=="toll_fam_mean" {
local xlab "xlab(-0.002(0.001)0.002)"
}	

if "`hhtype'"=="c" {
local fig a
}

if "`hhtype'"=="s" {
local fig b
}

coefplot ///
(`hhtype'_2008_2010_BEV_fam_yes, msym(O) mcol(navy) ciopts(lcol(navy))) (`hhtype'_2015_2017_BEV_fam_yes, msym(O) mcol(navy) ciopts(lcol(navy))) ///
(`hhtype'_2008_2010_car_fam_yes, msym(S) mcol(green) ciopts(lcol(green))) (`hhtype'_2015_2017_car_fam_yes, msym(S) mcol(green) ciopts(lcol(green))) ///
(`hhtype'_2008_2010_ICE_fam_count, msym(T) mcol(red) ciopts(lcol(red))) (`hhtype'_2015_2017_ICE_fam_count, msym(T) mcol(red) ciopts(lcol(red))) ///
(`hhtype'_2008_2010_car_fam_count, msym(D) mcol(brown) ciopts(lcol(brown))) (`hhtype'_2015_2017_car_fam_count, msym(D) mcol(brown) ciopts(lcol(brown))) ///
, keep(`treatvar'*)  aseq swapnames legend(off) `xlab' xline(0, lcol(gs6) lpattern(dash))  ///
headings(`hhtype'_2008_2010_BEV_fam_yes = "{bf:BEV (yes=1)}" `hhtype'_2008_2010_car_fam_yes = "{bf:Car (yes=1)}" `hhtype'_2008_2010_ICE_fam_count  = "{bf:ICEVs}" `hhtype'_2008_2010_car_fam_count = "{bf:Cars}"  ) ///
coeflabels(*2008_2010*="08-10"  *2015_2017*="15-17")  ///
graphregion(color(white)) plotregion(lcolor(black)) scale(1.3) 

graph export    "${figures}Figure4`fig'.pdf" ,  replace
graph export    "${figures}Figure4`fig'.png" ,  replace
}
}


